def sigmoid(z):
    return 1 / (1 + np.exp(-z))
#损失函数
def cost(theta, X, y, learningRate):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    first = np.multiply(-y, np.log(sigmoid(X * theta.T)))
    second = np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T)))
    reg = (learningRate / (2 * len(X))) * np.sum(np.power(theta[:,1:theta.shape[1]], 2))
    return np.sum(first - second) / len(X) + reg
#循环函数计算梯度

def gradient_with_loop(theta, X, y, learningRate):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    parameters = int(theta.ravel().shape[1])
    grad = np.zeros(parameters)
    
    error = sigmoid(X * theta.T) - y
    
    for i in range(parameters):
        term = np.multiply(error, X[:,i])
        
        if (i == 0):
            grad[i] = np.sum(term) / len(X)
        else:
            grad[i] = (np.sum(term) / len(X)) + ((learningRate / len(X)) * theta[:,i])
    
    return grad

#向量化梯度计算

def gradient(theta, X, y, learningRate):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    parameters = int(theta.ravel().shape[1])
    error = sigmoid(X * theta.T) - y
    
    # X转置后 400 * 5000 ，矩阵相加，对应元素相加，输出[400,]
    grad = ((X.T * error) / len(X)).T + ((learningRate / len(X)) * theta)
    
    # intercept gradient is not regularized
    grad[0, 0] = np.sum(np.multiply(error, X[:,0])) / len(X)
    
    return np.array(grad).ravel()
#需要计10个可能的类，并且由于逻辑回归只能一次在2个类之间进行分类，我们需要多类分类的策略。 其中具有k个不同类的标签就有k个分类器，每个分类器在“是类别 i”和“不是 i”之间决定，取值为 1  与 0 。 我们将把分类器训练包含在一个函数中，该函数计算10个分类器中的每个分类器的最终权重，并将权重返回为k ×（n + 1）数组，其中n是参数数量。

from scipy.optimize import minimize

def all_class_train(X, y, num_labels, learning_rate):
    rows = X.shape[0]
    params = X.shape[1]    
    # k* (n + 1) array for the parameters of each of the k classifiers   10 * 401
    all_theta = np.zeros((num_labels, params + 1))    
    # insert a column of ones at the beginning for the intercept term 在特征集的第一列插入1
    X = np.insert(X, 0, values=np.ones(rows), axis=1)    
    # labels are 1-indexed instead of 0-indexed  10个标签，分别计算10个逻辑回归模型
    for i in range(1, num_labels + 1):
        #401 个参数
        theta = np.zeros(params + 1)
        #转换成[1,0,0,10,0,0,0,0 .....]
        y_i = np.array([1 if label == i else 0 for label in y])
        #转置
        y_i = np.reshape(y_i, (rows, 1))        
        # minimize the objective function 
        fmin = minimize(fun=cost, x0=theta, args=(X, y_i, learning_rate), method='TNC', jac=gradient)
        #为当前分类标签，存储权重值
        all_theta[i-1,:] = fmin.x    
    return all_theta
rows = data['X'].shape[0]
params = data['X'].shape[1]

all_theta = np.zeros((10, params + 1))

X = np.insert(data['X'], 0, values=np.ones(rows), axis=1)

theta = np.zeros(params + 1)

y_0 = np.array([1 if label == 0 else 0 for label in data['y']])
y_0 = np.reshape(y_0, (rows, 1))

X.shape, y_0.shape, theta.shape, all_theta.shap
all_theta = all_class_train(data['X'], data['y'], 10, 1)
all_theta
#预测函数
def predict_all(X, all_theta):
    rows = X.shape[0]
    params = X.shape[1]
    num_labels = all_theta.shape[0]
    
    # same as before, insert ones to match the shape
    X = np.insert(X, 0, values=np.ones(rows), axis=1)
    
    # convert to matrices
    X = np.matrix(X)
    all_theta = np.matrix(all_theta)
    
    # compute the class probability for each class on each training instance  计算每个分类的概率
    h = sigmoid(X * all_theta.T)
    
    # create array of the index with the maximum probability  
    h_argmax = np.argmax(h, axis=1)
    
    # because our array was zero-indexed we need to add one for the true label prediction
    # 为了和 data['y']  中的数据比较
    h_argmax = h_argmax + 1
    
    return h_argmax
